##########################
# simulating IRRs (and risk differences)

## RUN LAST AFTER OTHER SCRIPTS -- note, takes a bit longer, maybe 3 mins?

# run to generate figures 2 and 3 -- simulations of 
# effects for subgroups

#this code generates an object called plot_ols and plot_k6
# which are th efigures for overall health and k6 respectively


#it also generates the estimates for the employment table
#########################




library(mvtnorm)
library(tidyverse)
library(cowplot)



#####################
# OLS Health Measure First
###################

#orig control group
survDesign <- svydesign(id      = ~0,
                        strata  = NULL,
                        weights = ~PERWEIGHT,
                        nest    = FALSE,
                        data    = sampledata2)

survDesign4.1<-subset(survDesign,!is.na(elig1))

(healthOLS <- svyglm(health ~ age_pol+ temp+SEX+REGION+YEAR+INTERVWMO+elig1*post, design=survDesign4.1))
summary(healthOLS, df.resid = degf(survDesign4.1))
confint(healthOLS)
nobs(healthOLS)

set.seed(02138)
sims<- rmvnorm(5000, mean = coef(healthOLS),
               sigma = vcov(healthOLS))


vals_data00 <- data.frame(model.matrix(healthOLS))
vals_data00<-vals_data00
vals_data00$post<-0
vals_data00$elig1<-0
vals_data00$elig1.post<-0
# vals_data$age_pol_cent<-i
vals_data00<-as.matrix(vals_data00)

X_beta <- t(vals_data00%*%t(sims))
expect_00 <-(X_beta)
expect_00_means<-rowMeans(expect_00)
mean(expect_00_means)


vals_data01 <- data.frame(model.matrix(healthOLS))
vals_data01<-vals_data01
vals_data01$post<-1
vals_data01$elig1<-0
vals_data01$elig1.post<-0
# vals_data$age_pol_cent<-i
vals_data01<-as.matrix(vals_data01)

X_beta <- t(vals_data01%*%t(sims))
expect_01 <-(X_beta)
expect_01_means<-rowMeans(expect_01)
mean(expect_01_means)

vals_data10 <- data.frame(model.matrix(healthOLS))
vals_data10<-vals_data10
vals_data10$post<-0
vals_data10$elig1<-1
vals_data10$elig1.post<-0
# vals_data$age_pol_cent<-i
vals_data10<-as.matrix(vals_data10)

X_beta <- t(vals_data10%*%t(sims))
expect_10 <-(X_beta)
expect_10_means<-rowMeans(expect_10)
mean(expect_10_means)


vals_data11 <- data.frame(model.matrix(healthOLS))
vals_data11<-vals_data11
vals_data11$post<-1
vals_data11$elig1<-1
vals_data11$elig1.post<-1
# vals_data$age_pol_cent<-i
vals_data11<-as.matrix(vals_data11)

X_beta <- t(vals_data11%*%t(sims))
expect_11 <-(X_beta)
expect_11_means<-rowMeans(expect_11)
mean(expect_11_means)


diffs<-(expect_11_means-expect_10_means)-(expect_01_means-expect_00_means)
mean(diffs)
quantile(diffs,c(0.025))
quantile(diffs,c(0.975))

group<-"cg1"
bindcg1<-data.frame(group,mean(diffs), quantile(diffs,c(0.025)),quantile(diffs,c(0.975)))



#group 2 -- met age criteria but not immig
survDesign <- svydesign(id      = ~0,
                        strata  = NULL,
                        weights = ~PERWEIGHT,
                        nest    = FALSE,
                        data    = sampledata2)

survDesign4.2<-subset(survDesign,!is.na(new_elig2))

(healthOLS <- svyglm(health ~ age_pol+ temp+SEX+REGION+YEAR+INTERVWMO+new_elig2*post, design=survDesign4.2))
summary(healthOLS, df.resid = degf(survDesign4.2))
confint(healthOLS)
nobs(healthOLS)


set.seed(02138)
sims<- rmvnorm(5000, mean = coef(healthOLS),
               sigma = vcov(healthOLS))


vals_data00 <- data.frame(model.matrix(healthOLS))
vals_data00<-vals_data00
vals_data00$post<-0
vals_data00$new_elig2<-0
vals_data00$new_elig2.post<-0
# vals_data$age_pol_cent<-i
vals_data00<-as.matrix(vals_data00)

X_beta <- t(vals_data00%*%t(sims))
expect_00 <-(X_beta)
expect_00_means<-rowMeans(expect_00)
mean(expect_00_means)


vals_data01 <- data.frame(model.matrix(healthOLS))
vals_data01<-vals_data01
vals_data01$post<-1
vals_data01$new_elig2<-0
vals_data01$new_elig2.post<-0
# vals_data$age_pol_cent<-i
vals_data01<-as.matrix(vals_data01)

X_beta <- t(vals_data01%*%t(sims))
expect_01 <-(X_beta)
expect_01_means<-rowMeans(expect_01)
mean(expect_01_means)

vals_data10 <- data.frame(model.matrix(healthOLS))
vals_data10<-vals_data10
vals_data10$post<-0
vals_data10$new_elig2<-1
vals_data10$new_elig2.post<-0
# vals_data$age_pol_cent<-i
vals_data10<-as.matrix(vals_data10)

X_beta <- t(vals_data10%*%t(sims))
expect_10 <-(X_beta)
expect_10_means<-rowMeans(expect_10)
mean(expect_10_means)


vals_data11 <- data.frame(model.matrix(healthOLS))
vals_data11<-vals_data11
vals_data11$post<-1
vals_data11$new_elig2<-1
vals_data11$new_elig2.post<-1
# vals_data$age_pol_cent<-i
vals_data11<-as.matrix(vals_data11)

X_beta <- t(vals_data11%*%t(sims))
expect_11 <-(X_beta)
expect_11_means<-rowMeans(expect_11)
mean(expect_11_means)


diffs<-(expect_11_means-expect_10_means)-(expect_01_means-expect_00_means)
mean(diffs)
quantile(diffs,c(0.025))
quantile(diffs,c(0.975))


group<-"cg2"
bindcg2<-data.frame(group,mean(diffs), quantile(diffs,c(0.025)),quantile(diffs,c(0.975)))




#group 3 - met immig criteria but not age
survDesign <- svydesign(id      = ~0,
                        strata  = NULL,
                        weights = ~PERWEIGHT,
                        nest    = FALSE,
                        data    = sampledata2)

survDesign4.3<-subset(survDesign,!is.na(new_elig3))

(healthOLS <- svyglm(health ~ new_elig3*post+age_pol+ temp+SEX+REGION+YEAR+INTERVWMO, design=survDesign4.3))
summary(healthOLS, df.resid = degf(survDesign4.3))
confint(healthOLS)
nobs(healthOLS)


set.seed(02138)
sims<- rmvnorm(5000, mean = coef(healthOLS),
               sigma = vcov(healthOLS))

vals_data00 <- data.frame(model.matrix(healthOLS))
vals_data00<-vals_data00
vals_data00$post<-0
vals_data00$new_elig3<-0
vals_data00$new_elig3.post<-0
# vals_data$age_pol_cent<-i
vals_data00<-as.matrix(vals_data00)

X_beta <- t(vals_data00%*%t(sims))
expect_00 <-(X_beta)
expect_00_means<-rowMeans(expect_00)
mean(expect_00_means)


vals_data01 <- data.frame(model.matrix(healthOLS))
vals_data01<-vals_data01
vals_data01$post<-1
vals_data01$new_elig3<-0
vals_data01$new_elig3.post<-0
# vals_data$age_pol_cent<-i
vals_data01<-as.matrix(vals_data01)

X_beta <- t(vals_data01%*%t(sims))
expect_01 <-(X_beta)
expect_01_means<-rowMeans(expect_01)
mean(expect_01_means)

vals_data10 <- data.frame(model.matrix(healthOLS))
vals_data10<-vals_data10
vals_data10$post<-0
vals_data10$new_elig3<-1
vals_data10$new_elig3.post<-0
# vals_data$age_pol_cent<-i
vals_data10<-as.matrix(vals_data10)

X_beta <- t(vals_data10%*%t(sims))
expect_10 <-(X_beta)
expect_10_means<-rowMeans(expect_10)
mean(expect_10_means)


vals_data11 <- data.frame(model.matrix(healthOLS))
vals_data11<-vals_data11
vals_data11$post<-1
vals_data11$new_elig3<-1
vals_data11$new_elig3.post<-1
# vals_data$age_pol_cent<-i
vals_data11<-as.matrix(vals_data11)

X_beta <- t(vals_data11%*%t(sims))
expect_11 <-(X_beta)
expect_11_means<-rowMeans(expect_11)
mean(expect_11_means)


diffs<-(expect_11_means-expect_10_means)-(expect_01_means-expect_00_means)
mean(diffs)
quantile(diffs,c(0.025))
quantile(diffs,c(0.975))

group<-"cg3"
bindcg3<-data.frame(group,mean(diffs), quantile(diffs,c(0.025)),quantile(diffs,c(0.975)))





#group 4 - met neither criteria
survDesign <- svydesign(id      = ~0,
                        strata  = NULL,
                        weights = ~PERWEIGHT,
                        nest    = FALSE,
                        data    = sampledata2)

survDesign4.4<-subset(survDesign,!is.na(new_elig4))

(healthOLS <- svyglm(health ~ age_pol+ temp+SEX+REGION+YEAR+INTERVWMO+new_elig4*post, design=survDesign4.4))
summary(healthOLS, df.resid = degf(survDesign4.4))
confint(healthOLS)
nobs(healthOLS)


set.seed(02138)
sims<- rmvnorm(5000, mean = coef(healthOLS),
               sigma = vcov(healthOLS))

vals_data00 <- data.frame(model.matrix(healthOLS))
vals_data00<-vals_data00
vals_data00$post<-0
vals_data00$new_elig4<-0
vals_data00$new_elig4.post<-0
# vals_data$age_pol_cent<-i
vals_data00<-as.matrix(vals_data00)

X_beta <- t(vals_data00%*%t(sims))
expect_00 <-(X_beta)
expect_00_means<-rowMeans(expect_00)
mean(expect_00_means)


vals_data01 <- data.frame(model.matrix(healthOLS))
vals_data01<-vals_data01
vals_data01$post<-1
vals_data01$new_elig4<-0
vals_data01$new_elig4.post<-0
# vals_data$age_pol_cent<-i
vals_data01<-as.matrix(vals_data01)

X_beta <- t(vals_data01%*%t(sims))
expect_01 <-(X_beta)
expect_01_means<-rowMeans(expect_01)
mean(expect_01_means)

vals_data10 <- data.frame(model.matrix(healthOLS))
vals_data10<-vals_data10
vals_data10$post<-0
vals_data10$new_elig4<-1
vals_data10$new_elig4.post<-0
# vals_data$age_pol_cent<-i
vals_data10<-as.matrix(vals_data10)

X_beta <- t(vals_data10%*%t(sims))
expect_10 <-(X_beta)
expect_10_means<-rowMeans(expect_10)
mean(expect_10_means)


vals_data11 <- data.frame(model.matrix(healthOLS))
vals_data11<-vals_data11
vals_data11$post<-1
vals_data11$new_elig4<-1
vals_data11$new_elig4.post<-1
# vals_data$age_pol_cent<-i
vals_data11<-as.matrix(vals_data11)

X_beta <- t(vals_data11%*%t(sims))
expect_11 <-(X_beta)
expect_11_means<-rowMeans(expect_11)
mean(expect_11_means)


diffs<-(expect_11_means-expect_10_means)-(expect_01_means-expect_00_means)
mean(diffs)
quantile(diffs,c(0.025))
quantile(diffs,c(0.975))




group<-"cg4"
bindcg4<-data.frame(group,mean(diffs), quantile(diffs,c(0.025)),quantile(diffs,c(0.975)))

labels<-c("1: R1. Didn't meet 1+ criteria","2: E1a. Met immig. age, but not policy age","3: E1b. Met policy age but not immig. age","4: E1c. Met neither criteria")
bind<-rbind(bindcg1,bindcg3,bindcg2,bindcg4)
bind<-cbind(labels,bind)

bind
plot_ols<-ggplot(data=bind,aes(c(1:4),y=mean.diffs.,group=labels,colour=labels))+ geom_hline(yintercept=0, size=1,linetype="dashed")+
  geom_errorbar(aes(ymin=quantile.diffs..c.0.025.., ymax=quantile.diffs..c.0.975..),  width=.15 , size=0.9) +geom_point(size=4)+
  theme(plot.title = element_text(hjust = 0.5))  + labs(x = 'Control group',  y = 'Effect (difference in difference in self-reported health)') +
  ggtitle("OLS-Self-reported Health") 

plot_ols




#####################################################################################################################

#poisson k6 

# sampledata2$age_pol<-as.factor(sampledata2$age_pol)
sampledata2$age_pol<-as.numeric(as.character((sampledata2$age_pol)))


#orig control group
survDesign2 <- svydesign(id      = ~0,
                         strata  = NULL,
                         weights = ~SAMPWEIGHT,
                         nest    = FALSE,
                         data    = sampledata2)

survDesign3.1<-subset(survDesign2,!is.na(mood_dist))

(k6POIS <- svyglm(k6_index ~ age_pol+ temp+SEX+REGION+YEAR+INTERVWMO+elig1*post, design=survDesign3.1, family=poisson))
summary(k6POIS, df.resid = degf(survDesign3.1))
confint(k6POIS)
nobs(k6POIS)





set.seed(02138)
sims<- rmvnorm(5000, mean = coef(k6POIS),
               sigma = vcov(k6POIS))

vals_data00 <- data.frame(model.matrix(k6POIS))
vals_data00<-vals_data00
vals_data00$post<-0
vals_data00$elig1<-0
vals_data00$elig1.post<-0
# vals_data$age_pol_cent<-i
vals_data00<-as.matrix(vals_data00)

X_beta <- t(vals_data00%*%t(sims))
expect_00 <-exp(X_beta)
# hist(expect_00)
# quantile(expect_00,c(0.99))
# 
# expect_00<-ifelse(expect_00>quantile(expect_00,c(0.99)),NA,expect_00)
expect_00_means<-rowMeans(expect_00,na.rm=TRUE)
mean(expect_00_means)


vals_data01 <- data.frame(model.matrix(k6POIS))
vals_data01<-vals_data01
vals_data01$post<-1
vals_data01$elig1<-0
vals_data01$elig1.post<-0
# vals_data$age_pol_cent<-i
vals_data01<-as.matrix(vals_data01)

X_beta <- t(vals_data01%*%t(sims))
expect_01 <-exp(X_beta)
# expect_01<-ifelse(expect_01>quantile(expect_01,c(0.99)),NA,expect_01)
expect_01_means<-rowMeans(expect_01,na.rm=TRUE)
mean(expect_01_means)






vals_data10 <- data.frame(model.matrix(k6POIS))
vals_data10<-vals_data10
vals_data10$post<-0
vals_data10$elig1<-1
vals_data10$elig1.post<-0
# vals_data$age_pol_cent<-i
vals_data10<-as.matrix(vals_data10)

X_beta <- t(vals_data10%*%t(sims))
expect_10 <-exp(X_beta)

# expect_10<-ifelse(expect_10>quantile(expect_10,c(0.99)),NA,expect_10)
expect_10_means<-rowMeans(expect_10,na.rm=TRUE)
mean(expect_10_means)






vals_data11 <- data.frame(model.matrix(k6POIS))
vals_data11<-vals_data11
vals_data11$post<-1
vals_data11$elig1<-1
vals_data11$elig1.post<-1
# vals_data$age_pol_cent<-i
vals_data11<-as.matrix(vals_data11)

X_beta <- t(vals_data11%*%t(sims))
expect_11 <-exp(X_beta)
# expect_11<-ifelse(expect_11>quantile(expect_11,c(0.99)),NA,expect_11)
expect_11_means<-rowMeans(expect_11,na.rm=TRUE)
mean(expect_11_means)



diffs<-(expect_11_means-expect_10_means)-(expect_01_means-expect_00_means)
mean(diffs)
quantile(diffs,c(0.025))
quantile(diffs,c(0.975))


IRRs<-(expect_11_means/expect_10_means)/(expect_01_means/expect_00_means)
mean(IRRs)
quantile(IRRs,c(0.025))
quantile(IRRs,c(0.975))


group<-"cg1"
bindcg1_k6<-data.frame(group,mean(diffs), quantile(diffs,c(0.025)),quantile(diffs,c(0.975)))




# s2<-subset(sampledata2,!is.na(mood_dist)& !is.na(new_elig2))
# table(s2$YEAR, s2$age_pol )


#group 2 -- met age criteria but not immig


# sampledata2$age_pol<-as.numeric(as.character(sampledata2$age_pol))
# sampledata2$age_pol<-as.factor(sampledata2$age_pol)

survDesign2 <- svydesign(id      = ~0,
                         strata  = NULL,
                         weights = ~SAMPWEIGHT,
                         nest    = FALSE,
                         data    = sampledata2)

survDesign3.2<-subset(survDesign2,!is.na(mood_dist) & !is.na(new_elig2))

(k6POIS <- svyglm(k6_index ~ age_pol+ temp+SEX+REGION+YEAR+INTERVWMO+new_elig2*post, design=survDesign3.2, family="poisson" ))
summary(k6POIS, df.resid = degf(survDesign3.2))
confint(k6POIS)
nobs(k6POIS)


set.seed(02138)
sims<- rmvnorm(5000, mean = coef(k6POIS),
               sigma = vcov(k6POIS))
vals_data00 <- data.frame(model.matrix(k6POIS))
colnames(vals_data00)
colnames(sims)

vals_data00<-vals_data00 
vals_data00$post<-0
vals_data00$new_elig2<-0
vals_data00$new_elig2.post<-0
# vals_data$age_pol_cent<-i
vals_data00<-as.matrix(vals_data00)

X_beta <- t(vals_data00%*%t(sims))
expect_00 <-exp(X_beta)
expect_00_means<-rowMeans(expect_00)
mean(expect_00_means)


vals_data01 <- data.frame(model.matrix(k6POIS))
vals_data01<-vals_data01
vals_data01$post<-1
vals_data01$new_elig2<-0
vals_data01$new_elig2.post<-0
# vals_data$age_pol_cent<-i
vals_data01<-as.matrix(vals_data01)

X_beta <- t(vals_data01%*%t(sims))
expect_01 <-exp(X_beta)
expect_01_means<-rowMeans(expect_01)
mean(expect_01_means)






vals_data10 <- data.frame(model.matrix(k6POIS))
vals_data10<-vals_data10 
vals_data10$post<-0
vals_data10$new_elig2<-1
vals_data10$new_elig2.post<-0
# vals_data$age_pol_cent<-i
vals_data10<-as.matrix(vals_data10)

X_beta <- t(vals_data10%*%t(sims))
expect_10 <-exp(X_beta)
expect_10_means<-rowMeans(expect_10)
mean(expect_10_means)




vals_data11 <- data.frame(model.matrix(k6POIS))
vals_data11<-vals_data11
vals_data11$post<-1
vals_data11$new_elig2<-1
vals_data11$new_elig2.post<-1
# vals_data$age_pol_cent<-i
vals_data11<-as.matrix(vals_data11)

X_beta <- t(vals_data11%*%t(sims))
expect_11 <-exp(X_beta)
expect_11_means<-rowMeans(expect_11)
mean(expect_11_means)



diffs<-(expect_11_means-expect_10_means)-(expect_01_means-expect_00_means)
mean(diffs)
quantile(diffs,c(0.025))
quantile(diffs,c(0.975))


IRRs<-(expect_11_means/expect_10_means)/(expect_01_means/expect_00_means)
mean(IRRs)
quantile(IRRs,c(0.025))
quantile(IRRs,c(0.975))


group<-"cg2"
bindcg2_k6<-data.frame(group,mean(diffs), quantile(diffs,c(0.025)),quantile(diffs,c(0.975)))









#group 3 - met immig criteria but not age
survDesign2 <- svydesign(id      = ~0,
                         strata  = NULL,
                         weights = ~SAMPWEIGHT,
                         nest    = FALSE,
                         data    = sampledata2)

survDesign3.3<-subset(survDesign2,!is.na(mood_dist) & !is.na(new_elig3))


(k6POIS <- svyglm(k6_index ~ new_elig3*post+age_pol+ temp+SEX+REGION+YEAR+INTERVWMO, design=survDesign3.3, family=poisson))
summary(k6POIS, df.resid = degf(survDesign3.3))
confint(k6POIS)
nobs(k6POIS)



set.seed(02138)
sims<- rmvnorm(5000, mean = coef(k6POIS),
               sigma = vcov(k6POIS))
vals_data00 <- data.frame(model.matrix(k6POIS))
vals_data00<-vals_data00
colnames(vals_data00)
vals_data00$post<-0
vals_data00$new_elig3<-0
vals_data00$new_elig3.post<-0
# vals_data$age_pol_cent<-i
vals_data00<-as.matrix(vals_data00)

X_beta <- t(vals_data00%*%t(sims))
expect_00 <-exp(X_beta)
expect_00_means<-rowMeans(expect_00)
mean(expect_00_means)


vals_data01 <- data.frame(model.matrix(k6POIS))
vals_data01<-vals_data01
vals_data01$post<-1
vals_data01$new_elig3<-0
vals_data01$new_elig3.post<-0
# vals_data$age_pol_cent<-i
vals_data01<-as.matrix(vals_data01)

X_beta <- t(vals_data01%*%t(sims))
expect_01 <-exp(X_beta)
expect_01_means<-rowMeans(expect_01)
mean(expect_01_means)






vals_data10 <- data.frame(model.matrix(k6POIS))
vals_data10<-vals_data10
vals_data10$post<-0
vals_data10$new_elig3<-1
vals_data10$new_elig3.post<-0
# vals_data$age_pol_cent<-i
vals_data10<-as.matrix(vals_data10)

X_beta <- t(vals_data10%*%t(sims))
expect_10 <-exp(X_beta)
expect_10_means<-rowMeans(expect_10)
mean(expect_10_means)




vals_data11 <- data.frame(model.matrix(k6POIS))
vals_data11<-vals_data11
vals_data11$post<-1
vals_data11$new_elig3<-1
vals_data11$new_elig3.post<-1
# vals_data$age_pol_cent<-i
vals_data11<-as.matrix(vals_data11)

X_beta <- t(vals_data11%*%t(sims))
expect_11 <-exp(X_beta)
expect_11_means<-rowMeans(expect_11)
mean(expect_11_means)

diffs<-(expect_11_means-expect_10_means)-(expect_01_means-expect_00_means)
mean(diffs)
quantile(diffs,c(0.025))
quantile(diffs,c(0.975))



IRRs<-(expect_11_means/expect_10_means)/(expect_01_means/expect_00_means)
mean(IRRs)
quantile(IRRs,c(0.025))
quantile(IRRs,c(0.975))


group<-"cg3"
bindcg3_k6<-data.frame(group,mean(diffs), quantile(diffs,c(0.025)),quantile(diffs,c(0.975)))



#group 4 - met neither criteria
survDesign2 <- svydesign(id      = ~0,
                         strata  = NULL,
                         weights = ~SAMPWEIGHT,
                         nest    = FALSE,
                         data    = sampledata2)

survDesign3.4<-subset(survDesign2,!is.na(mood_dist) & !is.na(new_elig4))


(k6POIS <- svyglm(k6_index ~ new_elig4*post+age_pol+ temp+SEX+REGION+YEAR+INTERVWMO, design=survDesign3.4, family=poisson))
summary(k6POIS, df.resid = degf(survDesign3.4))
confint(k6POIS)
nobs(k6POIS)


set.seed(02138)
sims<- rmvnorm(5000, mean = coef(k6POIS),
               sigma = vcov(k6POIS))
vals_data00 <- data.frame(model.matrix(k6POIS))
vals_data00<-vals_data00
colnames(vals_data00)
colnames(sims)

vals_data00$post<-0
vals_data00$new_elig4<-0
vals_data00$new_elig4.post<-0
# vals_data$age_pol_cent<-i
vals_data00<-as.matrix(vals_data00)

X_beta <- t(vals_data00%*%t(sims))
expect_00 <-exp(X_beta)
expect_00_means<-rowMeans(expect_00)
mean(expect_00_means)


vals_data01 <- data.frame(model.matrix(k6POIS))
vals_data01<-vals_data01
vals_data01$post<-1
vals_data01$new_elig4<-0
vals_data01$new_elig4.post<-0
# vals_data$age_pol_cent<-i
vals_data01<-as.matrix(vals_data01)

X_beta <- t(vals_data01%*%t(sims))
expect_01 <-exp(X_beta)
expect_01_means<-rowMeans(expect_01)
mean(expect_01_means)






vals_data10 <- data.frame(model.matrix(k6POIS))
vals_data10<-vals_data10
vals_data10$post<-0
vals_data10$new_elig4<-1
vals_data10$new_elig4.post<-0
# vals_data$age_pol_cent<-i
vals_data10<-as.matrix(vals_data10)

X_beta <- t(vals_data10%*%t(sims))
expect_10 <-exp(X_beta)
expect_10_means<-rowMeans(expect_10)
mean(expect_10_means)




vals_data11 <- data.frame(model.matrix(k6POIS))
vals_data11<-vals_data11
vals_data11$post<-1
vals_data11$new_elig4<-1
vals_data11$new_elig4.post<-1
# vals_data$age_pol_cent<-i
vals_data11<-as.matrix(vals_data11)

X_beta <- t(vals_data11%*%t(sims))
expect_11 <-exp(X_beta)
expect_11_means<-rowMeans(expect_11)
mean(expect_11_means)



diffs<-(expect_11_means-expect_10_means)-(expect_01_means-expect_00_means)
mean(diffs)
quantile(diffs,c(0.025))
quantile(diffs,c(0.975))



IRRs<-(expect_11_means/expect_10_means)/(expect_01_means/expect_00_means)
mean(IRRs)
quantile(IRRs,c(0.025))
quantile(IRRs,c(0.975))



group<-"cg4"
bindcg4_k6<-data.frame(group,mean(diffs), quantile(diffs,c(0.025)),quantile(diffs,c(0.975)))



labels<-c("1: R1. Didn't meet 1+ criteria","2: E1a. Met immig. age, but not policy age","3: E1b. Met policy age but not immig. age","4: E1c. Met neither criteria")
bind_k6<-rbind(bindcg1_k6,bindcg3_k6,bindcg2_k6,bindcg4_k6)
bind_k6<-cbind(labels,bind_k6)


plot_k6<-ggplot(data=bind_k6,aes(c(1:4),y=mean.diffs.,group=labels,colour=labels))+ geom_hline(yintercept=0, size=1,linetype="dashed")+
  geom_errorbar(aes(ymin=quantile.diffs..c.0.025.., ymax=quantile.diffs..c.0.975..),  width=.15 , size=0.9) +geom_point(size=4)+
  theme(plot.title = element_text(hjust = 0.5))  + labs(x = 'Control group',  y = 'Effect (difference in difference in K6 score)') +
  ggtitle("Poisson-K6 Score")

plot_k6



###############################################################################
# employment

survDesign2 <- svydesign(id      = ~0,
                         strata  = NULL,
                         weights = ~SAMPWEIGHT,
                         nest    = FALSE,
                         data    = sampledata2)

survDesign3.1<-subset(survDesign2,!is.na(emp_year))

employedLOGIS <- svyglm(emp_year ~ AGE+temp+SEX+REGION+YEAR+INTERVWMO+elig1*post, design=survDesign3.1, family=binomial (link="logit"))
summary(employedLOGIS, df.resid = degf(survDesign3.1))
confint(employedLOGIS)
nobs(employedLOGIS)




set.seed(02138)
sims<- rmvnorm(5000, mean = coef(employedLOGIS),
               sigma = vcov(employedLOGIS))

vals_data00 <- data.frame(model.matrix(employedLOGIS))
vals_data00<-vals_data00
vals_data00$post<-0
vals_data00$elig1<-0
vals_data00$elig1.post<-0
# vals_data$age_pol_cent<-i
vals_data00<-as.matrix(vals_data00)

X_beta <- t(vals_data00%*%t(sims))
expect_00 <-exp(X_beta)/(1+exp(X_beta))
# hist(expect_00)
# quantile(expect_00,c(0.99))
# 
# expect_00<-ifelse(expect_00>quantile(expect_00,c(0.99)),NA,expect_00)
expect_00_means<-rowMeans(expect_00,na.rm=TRUE)
mean(expect_00_means)


vals_data01 <- data.frame(model.matrix(employedLOGIS))
vals_data01<-vals_data01
vals_data01$post<-1
vals_data01$elig1<-0
vals_data01$elig1.post<-0
# vals_data$age_pol_cent<-i
vals_data01<-as.matrix(vals_data01)

X_beta <- t(vals_data01%*%t(sims))
expect_01 <-exp(X_beta)/(1+exp(X_beta))
# expect_01<-ifelse(expect_01>quantile(expect_01,c(0.99)),NA,expect_01)
expect_01_means<-rowMeans(expect_01,na.rm=TRUE)
mean(expect_01_means)






vals_data10 <- data.frame(model.matrix(employedLOGIS))
vals_data10<-vals_data10
vals_data10$post<-0
vals_data10$elig1<-1
vals_data10$elig1.post<-0
# vals_data$age_pol_cent<-i
vals_data10<-as.matrix(vals_data10)

X_beta <- t(vals_data10%*%t(sims))
expect_10 <-exp(X_beta)/(1+exp(X_beta))

# expect_10<-ifelse(expect_10>quantile(expect_10,c(0.99)),NA,expect_10)
expect_10_means<-rowMeans(expect_10,na.rm=TRUE)
mean(expect_10_means)






vals_data11 <- data.frame(model.matrix(employedLOGIS))
vals_data11<-vals_data11
vals_data11$post<-1
vals_data11$elig1<-1
vals_data11$elig1.post<-1
# vals_data$age_pol_cent<-i
vals_data11<-as.matrix(vals_data11)

X_beta <- t(vals_data11%*%t(sims))
expect_11 <-exp(X_beta)/(1+exp(X_beta))
# expect_11<-ifelse(expect_11>quantile(expect_11,c(0.99)),NA,expect_11)
expect_11_means<-rowMeans(expect_11,na.rm=TRUE)
mean(expect_11_means)



diffs<-(expect_11_means-expect_10_means)-(expect_01_means-expect_00_means)
mean(diffs)
quantile(diffs,c(0.025))
quantile(diffs,c(0.975))


# ORs<-(expect_11_means/expect_10_means)/(expect_01_means/expect_00_means)
# mean(IRRs)
# quantile(IRRs,c(0.025))
# quantile(IRRs,c(0.975))


group<-"cg1"
bindcg1_em<-data.frame(group,mean(diffs), quantile(diffs,c(0.025)),quantile(diffs,c(0.975)))




# s2<-subset(sampledata2,!is.na(poor_health)& !is.na(new_elig2))
# table(s2$YEAR, s2$age_pol )





#group 2 -- met age criteria but not immig


# sampledata2$age_pol<-as.numeric(as.character(sampledata2$age_pol))
# sampledata2$age_pol<-as.factor(sampledata2$age_pol)

survDesign3.2<-subset(survDesign2,!is.na(emp_year) & !is.na(new_elig2))

employedLOGIS <- svyglm(emp_year ~ AGE+temp+SEX+REGION+YEAR+INTERVWMO+new_elig2*post, design=survDesign3.2, family=binomial (link="logit"))
summary(employedLOGIS, df.resid = degf(survDesign3.2))
confint(employedLOGIS)
nobs(employedLOGIS) 



set.seed(02138)
sims<- rmvnorm(5000, mean = coef(employedLOGIS),
               sigma = vcov(employedLOGIS))

vals_data00 <- data.frame(model.matrix(employedLOGIS))
vals_data00<-vals_data00
vals_data00$post<-0
vals_data00$new_elig2<-0
vals_data00$new_elig2.post<-0
# vals_data$age_pol_cent<-i
vals_data00<-as.matrix(vals_data00)

X_beta <- t(vals_data00%*%t(sims))
expect_00 <-exp(X_beta)/(1+exp(X_beta))
# hist(expect_00)
# quantile(expect_00,c(0.99))
# 
# expect_00<-ifelse(expect_00>quantile(expect_00,c(0.99)),NA,expect_00)
expect_00_means<-rowMeans(expect_00,na.rm=TRUE)
mean(expect_00_means)


vals_data01 <- data.frame(model.matrix(employedLOGIS))
vals_data01<-vals_data01
vals_data01$post<-1
vals_data01$new_elig2<-0
vals_data01$new_elig2.post<-0
# vals_data$age_pol_cent<-i
vals_data01<-as.matrix(vals_data01)

X_beta <- t(vals_data01%*%t(sims))
expect_01 <-exp(X_beta)/(1+exp(X_beta))
# expect_01<-ifelse(expect_01>quantile(expect_01,c(0.99)),NA,expect_01)
expect_01_means<-rowMeans(expect_01,na.rm=TRUE)
mean(expect_01_means)






vals_data10 <- data.frame(model.matrix(employedLOGIS))
vals_data10<-vals_data10
vals_data10$post<-0
vals_data10$new_elig2<-1
vals_data10$new_elig2.post<-0
# vals_data$age_pol_cent<-i
vals_data10<-as.matrix(vals_data10)

X_beta <- t(vals_data10%*%t(sims))
expect_10 <-exp(X_beta)/(1+exp(X_beta))

# expect_10<-ifelse(expect_10>quantile(expect_10,c(0.99)),NA,expect_10)
expect_10_means<-rowMeans(expect_10,na.rm=TRUE)
mean(expect_10_means)






vals_data11 <- data.frame(model.matrix(employedLOGIS))
vals_data11<-vals_data11
vals_data11$post<-1
vals_data11$new_elig2<-1
vals_data11$new_elig2.post<-1
# vals_data$age_pol_cent<-i
vals_data11<-as.matrix(vals_data11)

X_beta <- t(vals_data11%*%t(sims))
expect_11 <-exp(X_beta)/(1+exp(X_beta))
# expect_11<-ifelse(expect_11>quantile(expect_11,c(0.99)),NA,expect_11)
expect_11_means<-rowMeans(expect_11,na.rm=TRUE)
mean(expect_11_means)



diffs<-(expect_11_means-expect_10_means)-(expect_01_means-expect_00_means)
mean(diffs)
quantile(diffs,c(0.025))
quantile(diffs,c(0.975))


# ORs<-(expect_11_means/expect_10_means)/(expect_01_means/expect_00_means)
# mean(IRRs)
# quantile(IRRs,c(0.025))
# quantile(IRRs,c(0.975))


group<-"cg2"
bindcg2_em<-data.frame(group,mean(diffs), quantile(diffs,c(0.025)),quantile(diffs,c(0.975)))









#group 3 - met immig criteria but not age
survDesign3.3<-subset(survDesign2,!is.na(emp_year) & !is.na(new_elig3))

employedLOGIS <- svyglm(emp_year ~ AGE+temp+SEX+REGION+YEAR+INTERVWMO+new_elig3*post, design=survDesign3.3, family=binomial (link="logit"))
summary(employedLOGIS, df.resid = degf(survDesign3.3))
confint(employedLOGIS)
nobs(employedLOGIS)






set.seed(02138)
sims<- rmvnorm(5000, mean = coef(employedLOGIS),
               sigma = vcov(employedLOGIS))

vals_data00 <- data.frame(model.matrix(employedLOGIS))
vals_data00<-vals_data00
vals_data00$post<-0
vals_data00$new_elig3<-0
vals_data00$new_elig3.post<-0
# vals_data$age_pol_cent<-i
vals_data00<-as.matrix(vals_data00)

X_beta <- t(vals_data00%*%t(sims))
expect_00 <-exp(X_beta)/(1+exp(X_beta))
# hist(expect_00)
# quantile(expect_00,c(0.99))
# 
# expect_00<-ifelse(expect_00>quantile(expect_00,c(0.99)),NA,expect_00)
expect_00_means<-rowMeans(expect_00,na.rm=TRUE)
mean(expect_00_means)


vals_data01 <- data.frame(model.matrix(employedLOGIS))
vals_data01<-vals_data01
vals_data01$post<-1
vals_data01$new_elig3<-0
vals_data01$new_elig3.post<-0
# vals_data$age_pol_cent<-i
vals_data01<-as.matrix(vals_data01)

X_beta <- t(vals_data01%*%t(sims))
expect_01 <-exp(X_beta)/(1+exp(X_beta))
# expect_01<-ifelse(expect_01>quantile(expect_01,c(0.99)),NA,expect_01)
expect_01_means<-rowMeans(expect_01,na.rm=TRUE)
mean(expect_01_means)






vals_data10 <- data.frame(model.matrix(employedLOGIS))
vals_data10<-vals_data10
vals_data10$post<-0
vals_data10$new_elig3<-1
vals_data10$new_elig3.post<-0
# vals_data$age_pol_cent<-i
vals_data10<-as.matrix(vals_data10)

X_beta <- t(vals_data10%*%t(sims))
expect_10 <-exp(X_beta)/(1+exp(X_beta))

# expect_10<-ifelse(expect_10>quantile(expect_10,c(0.99)),NA,expect_10)
expect_10_means<-rowMeans(expect_10,na.rm=TRUE)
mean(expect_10_means)






vals_data11 <- data.frame(model.matrix(employedLOGIS))
vals_data11<-vals_data11
vals_data11$post<-1
vals_data11$new_elig3<-1
vals_data11$new_elig3.post<-1
# vals_data$age_pol_cent<-i
vals_data11<-as.matrix(vals_data11)

X_beta <- t(vals_data11%*%t(sims))
expect_11 <-exp(X_beta)/(1+exp(X_beta))
# expect_11<-ifelse(expect_11>quantile(expect_11,c(0.99)),NA,expect_11)
expect_11_means<-rowMeans(expect_11,na.rm=TRUE)
mean(expect_11_means)



diffs<-(expect_11_means-expect_10_means)-(expect_01_means-expect_00_means)
mean(diffs)
quantile(diffs,c(0.025))
quantile(diffs,c(0.975))


# ORs<-(expect_11_means/expect_10_means)/(expect_01_means/expect_00_means)
# mean(IRRs)
# quantile(IRRs,c(0.025))
# quantile(IRRs,c(0.975))


group<-"cg3"
bindcg3_em<-data.frame(group,mean(diffs), quantile(diffs,c(0.025)),quantile(diffs,c(0.975)))




#group 4 - met neither criteria

survDesign3.4<-subset(survDesign2,!is.na(emp_year) & !is.na(new_elig4))

employedLOGIS <- svyglm(emp_year ~ AGE+temp+SEX+REGION+YEAR+INTERVWMO+new_elig4*post, design=survDesign3.4, family=binomial (link="logit"))
summary(employedLOGIS, df.resid = degf(survDesign3.4))
confint(employedLOGIS)
nobs(employedLOGIS)

set.seed(02138)
sims<- rmvnorm(5000, mean = coef(employedLOGIS),
               sigma = vcov(employedLOGIS))

vals_data00 <- data.frame(model.matrix(employedLOGIS))
vals_data00<-vals_data00
vals_data00$post<-0
vals_data00$new_elig4<-0
vals_data00$new_elig4.post<-0
# vals_data$age_pol_cent<-i
vals_data00<-as.matrix(vals_data00)

X_beta <- t(vals_data00%*%t(sims))
expect_00 <-exp(X_beta)/(1+exp(X_beta))
# hist(expect_00)
# quantile(expect_00,c(0.99))
# 
# expect_00<-ifelse(expect_00>quantile(expect_00,c(0.99)),NA,expect_00)
expect_00_means<-rowMeans(expect_00,na.rm=TRUE)
mean(expect_00_means)


vals_data01 <- data.frame(model.matrix(employedLOGIS))
vals_data01<-vals_data01
vals_data01$post<-1
vals_data01$new_elig4<-0
vals_data01$new_elig4.post<-0
# vals_data$age_pol_cent<-i
vals_data01<-as.matrix(vals_data01)

X_beta <- t(vals_data01%*%t(sims))
expect_01 <-exp(X_beta)/(1+exp(X_beta))
# expect_01<-ifelse(expect_01>quantile(expect_01,c(0.99)),NA,expect_01)
expect_01_means<-rowMeans(expect_01,na.rm=TRUE)
mean(expect_01_means)






vals_data10 <- data.frame(model.matrix(employedLOGIS))
vals_data10<-vals_data10
vals_data10$post<-0
vals_data10$new_elig4<-1
vals_data10$new_elig4.post<-0
# vals_data$age_pol_cent<-i
vals_data10<-as.matrix(vals_data10)

X_beta <- t(vals_data10%*%t(sims))
expect_10 <-exp(X_beta)/(1+exp(X_beta))

# expect_10<-ifelse(expect_10>quantile(expect_10,c(0.99)),NA,expect_10)
expect_10_means<-rowMeans(expect_10,na.rm=TRUE)
mean(expect_10_means)






vals_data11 <- data.frame(model.matrix(employedLOGIS))
vals_data11<-vals_data11
vals_data11$post<-1
vals_data11$new_elig4<-1
vals_data11$new_elig4.post<-1
# vals_data$age_pol_cent<-i
vals_data11<-as.matrix(vals_data11)

X_beta <- t(vals_data11%*%t(sims))
expect_11 <-exp(X_beta)/(1+exp(X_beta))
# expect_11<-ifelse(expect_11>quantile(expect_11,c(0.99)),NA,expect_11)
expect_11_means<-rowMeans(expect_11,na.rm=TRUE)
mean(expect_11_means)



diffs<-(expect_11_means-expect_10_means)-(expect_01_means-expect_00_means)
mean(diffs)
quantile(diffs,c(0.025))
quantile(diffs,c(0.975))


# ORs<-(expect_11_means/expect_10_means)/(expect_01_means/expect_00_means)
# mean(IRRs)
# quantile(IRRs,c(0.025))
# quantile(IRRs,c(0.975))


group<-"cg4"
bindcg4_em<-data.frame(group,mean(diffs), quantile(diffs,c(0.025)),quantile(diffs,c(0.975)))





labels<-c("1: R1. Didn't meet 1+ criteria","2: E1a. Met immig. age, but not policy age","3: E1b. Met policy age but not immig. age","4: E1c. Met neither criteria")
bind_em<-rbind(bindcg1_em,bindcg3_em,bindcg2_em,bindcg4_em)
bind_em<-cbind(labels,bind_em)
bind_em


plot_em<-ggplot(data=bind_em,aes(c(1:4),y=mean.diffs.,group=labels,colour=labels))+ geom_hline(yintercept=0, size=1,linetype="dashed")+
  geom_errorbar(aes(ymin=quantile.diffs..c.0.025.., ymax=quantile.diffs..c.0.975..),  width=.15 , size=0.9) +geom_point(size=4)+
  theme(plot.title = element_text(hjust = 0.5))  + labs(x = 'Control group',  y = 'Effect (difference in difference in P(employed))') +
  ggtitle("P(employed)")

plot_em



